Active Brownian Motion in Threshold Distribution of a Coulomb Blockade Model 
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Randomly-distributed offset charges affect the nonlinear current-voltage property via the fluctuation of the 
threshold voltage of Coulomb blockade arrays. We analytically derive the distribution of the threshold voltage 
for a model of one-dimensional locally-coupled Coulomb blockade arrays, and propose a general relationship 
between conductance and the distribution. In addition, we show the distribution for a long array is equivalent to 
the distribution of the number of upward steps for aligned objects of different height. The distribution satisfies 
a novel Fokker-Planck equation corresponding to active Brownian motion. The feature of the distribution is 
clarified by comparing it with the Wigner and Ornstein-Uhlenbeck processes. It is not restricted to the Coulomb 
blockade model, but instructive in statistical physics generally. 
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Introduction. — Nonlinear phenomena and threshold behav- 
iors are observed in many disordered systems QjJ] . A Coulomb 
blockade (CB) [2] is one such example for which charac- 
teristically nonlinear current-voltage (I-V) behavior occurs 
above a threshold voltage V±. Specifically, CB is the in- 
creased resistance at low bias voltage of an electronic device 
having a low-capacitance tunnel junction, the thin insulating 
barrier that lies between two electrodes across which electrons 
tunnel quantum mechanically. Owing to CB, the conductance 
of the device is not constant at low voltage, and no current 
flows below VJh. 

Studies have explicitly considered types of disorder and 
clarified that disorder affects transport phenomena 
Middleton and Wingreen (MW) considered the charge dis- 
order that originates from impurities of a substrate yfl . The 
threshold voltage is sensitive to this charge disorder. The dis- 
tribution of V± has never been derived, although MW have 
discussed the mean value and variance |3,|g]. 

In this Letter, we focus on the threshold distribution (TD) 
as it leads to understanding the nonlinearity in I-V response; 
we show that the conductance is represented by the cumula- 
tive distribution of V±. We find an analytic expression for 
the TD for a one-dimensional (ID) locally coupled CB array. 
In addition, we reveal that the TD in the long-array limit is 
equivalent to the distribution for the number of upward steps 
for aligned objects of different height. The distribution satis- 
fies a novel Fokker-Planck equation corresponding to active 
Brownian motion |8|]; i.e., overdamped motion of a Brownian 
particle in a harmonic potential that spreads with time. This 
characteristic of the distribution is quite instructive in the field 
of statistical physics. 

Model. — We employ the model proposed by MW 13], in 
which there are TV aligned Coulomb islands, constituting the 
minimum units of charge storage (Fig. We consider that 
the gate capacitance C 9 is much greater than the island-island 
and island-electrode capacitances C. In general, interactions 
such as electron-electron and spin-coupling play an impor- 
tant role in evolving the nonlinear I-V behavior |a,|2D. How- 
ever, such interactions are not dominant if C/C 9 <C 1 corre- 
sponding to the so-called locally-coupled CB. Compared with 
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FIG. 1. A configuration of a ID array with N = 8 (upper) and a 
distribution of the offset charge q\ (bottom). In the upper figure, a cir- 
cle represents a Coulomb island and the connecting line a tunneling 
junction with capacitance C. The array is sandwiched between posi- 
tive (+) and negative (— ) electrodes. Each island connects to a gate 
electrode, omitted in the figure, with capacitance C 9 . In the bottom 
figure, an arrow indicates an upward step, which means q\ < gi+i 
(1 G {1, . . . , N — 1}). In this distribution, there are three upward 
steps and four offset charges less than q± . 



several theoretical approaches such as the density functional 
theory @, [Hi and the random matrix theory fllll , the model 
we employ is classical and the simplest demonstrating CB. 
Many experimental features though can be explained by this 
model T2UT5I1 . and theoretical work is still continuing even 
7.ll6ll. More significantly, some results obtained in 
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this Letter are not solely restricted to this model. 

The voltages of the negative and gate electrodes are set to 
zero, and the bias voltage is thus equivalent to the voltage 
<£>+ of the positive electrode. Let Qi denote the charge of 
the i-th island; i G {1, . . . , N}. The charge is represented 
as Qi = n e e + q x , where n e denotes an integer, e the ele- 
mentary charge, and qi the offset charge arise from an impu- 
rity. The offset charges are given by uniform random num- 
bers in [— e/2, e/2], and remain constant over time. The off- 
set charges just indicate the non-integral part of each charge, 
i.e., the uniform distribution for q x is equivalent to arbitrary 
distributions of offset charges. 

The total energy E of the system is written as 
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where Q+ denotes the charge of the positive electrode. My 
denotes the capacitance matrix; for ID simple arrays, My = 
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C g + 2C for i = j, Mj = -C for |i - j| = 1, and Mj = 
otherwise. The system evolves such that E decreases. To take 
the most probable path of evolution, we transfer an electron to 
another island and calculate the energy change AE^^y for all 
possible tunneling paths, where {i', j'} G {1, . . . , N, +, — }. 
In simulations (e.g., USES), each tunneling time, which 
is proportional to the change in energy for T = iQ, is 
calculated, and the shortest tunneling time is thus employed 
for the time evolution increments. In the rest of the paper, we 
work in dimensionless units whereby the charge is scaled by 
e, the voltage by e/C g , and the energy by e 2 /C g . 

Vth as a function of q\ . — As a simple example, let us con- 
sider an array with N = 2 and describe V± as a function 
of offset charges. There are six possible paths; however, it 

is sufficient to consider A£?i^ + , AE 2 ^i, and AE >2 for 

<£ + > 0. Note that the paths in the reverse direction should be 
considered when <£> + < 0. In the limit C/C g — >• 0, 



AE^+ < 
AE 2 ^i < ^ 
AE_^2 < ^ Q 2 > 1/2 



>Qi + l/2, 
Qi-Q2> l, 



(2a) 
(2b) 
(2c) 



If all energy changes are greater than zero, no electrons get 
transferred; i.e., blockading occurs. Equation (l2bl) suggests 
that it is effective to separately consider charge-offset condi- 
tions qi > q2 (no upward steps) and q\ < q2 (an upward 
step). As <1> + increases quasi- statically, under the former con- 
dition, Eq. (l2ab is satisfied above <£ + = q\ + 1/2, and an elec- 
tron then is transferred from island 1 to the positive electrode. 
Thus, Eq. (l2bl) and subsequently Eq. (l2cl are satisfied. After- 
ward, Eq. (l2at is again satisfied. This cycle consequently gets 
repeated; i.e., the current flows between the positive and neg- 
ative electrodes in a steady state above $ + > Vth = gi + 1/2. 
In contrast, in the latter case, even if Eq. (l2at is satisfied and 
an electron moves from island 1 to the positive electrode, 
AE 2 ^\ remains greater than zero because q\ < q 2 . For 
AEi^ + and AE 2 ^\ to be less than zero, <1> + has to be in- 
creased to qi + 3/2, and a steady-state current then flows; i.e., 
the voltage threshold is V± = gi + 3/2. The above argument 
holds, without loss of generality, to arbitrary N; i.e., 

V th (q u n) = qi +n-l/2 (-1/2 < q x < 1/2) (3a) 
& qi(V th ,n) = V th -n + l/2 (n-l<V th <n) (3b) 

where n—1 indicates the number of upward steps; 1 < n < 
N. The threshold depends only on q\ and n; i.e., the mag- 
nitudes of the offset charges between neighboring islands is 
renormalized to n. 

Threshold distribution. — Equation (l3at suggests that the 
charge-offset analysis based on q\ is appropriate. In addi- 
tion, Eq. (l3bl) suggests that the region < V± < N should 
be divided into N equally- spaced segments. Thus, the n-th 
segmented TD for the -island array is expressed as 
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U N (n\k)U N (k) (n-l<V th <n), (4) 



where Ujsf(n\k) denotes the conditional probability that there 
are n — 1 upward steps if there are k offset charges less than 
<7h' • Note that U]\r(n\k) does not depend on V±. Here, since q\ 
is the basis for analyzing the offset charges, we should select 
f/ = 1. Il v(fc) denotes the probability that there are k offset 
charges less than q\ , and is expressed as 
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where pc and pl are the probabilities of q^ > q\ and q^ < q\ , 
respectively, and h G {2, . . . , AT}. Note that po = 1/2 — q\ 
and p L = 1/2 + qi. 

One can obtain U 2 (1\0) = U 2 (2\l) = and U 2 (l\l) = 
U 2 (2\0) = 1, and then, 



P 2 (1) (V th ) = V th , P^(V th ) = 2-V t] 



(2)/ 



(6) 



Using the same procedure, we obtain the entire TD Pn(V±) 
for arbitrary N as the joining of the segmented TDs pj^ (V±) 
II19I1 . As shown in Fig. [3 simulation results are correctly de- 
scribed without fitting parameters. It is clear that, for arbi- 
trary N, each segmented TD is represented as an (N — 1)- 
degree polynomial of Vth because of the term PL k PG N ~ 1 ~ k • 
For small N (in particular, N = 2 in Fig.0, the distributions 
have strange shape which might be a consequence of model- 
dependent behavior. In more realistic cases, other physical 
effects such as electrode shape should be taken into account. 

Distribution of upward steps. — The conditional probabil- 
ity Ujsr(n\k) determines the TD for arbitrary N. However, in 
practice, it is difficult to obtain UN(n\k) for large N. To in- 
vestigate the TD for large N, we focus on the intersections of 
the segmented TDs. In particular, we focus on the right edge 
of each segment; i.e., Vth = n- Since (pl,Pg) = (1, 0) at the 
right edge, Eq. dU reduces to 
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U N (n\N-l)=:Y(n,N) (atVfc = n). (7) 



Therefore, our problem results in obtaining Y(n, N) that in- 
dicates the probability in the case of n — 1 upward steps for 
N — 1 aligned objects (i.e., q 2 , . . . , qw) of different height. 
Since none of the specific features of the model are used, the 
discussion in the rest of this section is not limited to CB but 
has applicability to statistical physics generally. 

We consider the probability that the number of upward 
steps for N + 1 different heights is the same as that for N 
different heights. According to Fig. the probability is ex- 
pressed by (k + 1) /N, where the brackets (•) indicate the av- 
erage for 
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^tvIi(^) denotes the probability that there are k offset 
charges less than gy if there are n — 1 upward steps in A^ — 1 
offset charges, where the basis for analyzing offset charges is 
qw, i.e., f/ = N. Although a mathematical proof has yet to be 
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FIG. 2. (Color online) Plot of TDs of Vth for N = 2, 3, 4, and 5 (from right to left). A filled circle represents a simulation result and a colored 
line the segmented TD P^(Vth) obtained analytically for n — 1 (blue), n = 2 (green), n — 3 (brown), n — 4 (red), and n — 5 (purple). 

(Vth) for iV = 2 is expressed by Eq. J6]) and expressions for N = 3, 4, and 5 are given in the supplement fl9h . The inset is a close-up 
(semi-log plot) of the first segment for N = 5. The simulation used 10 6 different initial distributions of the offset charges. 
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TABLE I. Comparison of the variance corresponding to Eq. l fT2t . 
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FIG. 3. An example of the distribution, where there are five upward 
steps (n = 6) and five offset charges less than qN (k = 5). The 
arrows at right indicate the possible location of gjv+i. There are N 
arrows in total, with k + 1 arrows belonging to (A). If an arrow is 
chosen from (A), no increase in upward steps occurs. 

given, the probability (k + 1) /N is expected to be n/N lEoll . 
This expectation is understandable qualitatively as follows. If 
there are already many upward steps (i.e., large n), then q^ 
tends to be greater than other offset charges (#2, • • • , Qn-i)- 
Thus, the probability tends to increase with increasing n. With 
this expectation, the recurrence formula for Y(n,N) is ob- 
tained as 

n at _ ( n _ i ) 
F(n,7V+l) = -F(n,7V) + y - J -Y{n-1,N). (9) 

For later discussion, we introduce both a Active field x = 
n — N/2 and time t = N. Note that x and t do not indicate 
electron motions, but are just changes in variables [21]. By 
defining Z(x, t) := Y(x + AT/2, TV), Eq. ® reduces to 

Z(M+1) = Q + ^) Z(x+,t)+Q - ^jz(x-,t) , 

(10) 

where x± = x ± 1/2. In the continuous limit, a partial differ- 
ential equation is obtained 

dZ(x,t) = d_rx 1 &Z{x,t) 

dt dxlt [ '>\ + dx* ' ( } 

for which D = 1/8 describes locally-coupled CB. The first 
term of r.h.s. depends explicitly on time, so that the equation is 



classified as related to a time-dependent Ornstein-Uhlenbeck 
(OU) process ll22ll . The differential equation is equivalent to 
the Fokker-Planck equation corresponding to active Brownian 
motion [8]; i.e., overdamped motion of a Brownian particle in 
a harmonic potential <fr(X, t) = X 2 /2t, represented as 

-^ = -^<l>{X,t) + V2Di{t), (12) 

where X = X(t) denotes the position of the Brownian 
particle, and denotes a fluctuating term that satisfies 
= 5 (t ~ O with delta Unction S(-). This novel 
relationship between the distribution of the upward steps and 
the active Brownian motion is analogous to that between the 
binomial coefficient and Brownian motion. 

It can be shown that the distribution Z(x,t) is Gaussian 
with variance 3Dt/2 under the limit t —> 00 112311 . In that 
limit, although the variance of the OU process (i.e., t) = 
X 2 /2 in Eq. (fT2b ) is a constant D, that of the above time- 
dependent OU process is proportional to t (Tabled]). This is 
qualitatively the same as the Wiener process (i.e., <fr(X, t) = 
in Eq. (fT2b): however, the variance of Z(x, t) is smaller than 
that of the Wiener process of 2Dt. The presence of the poten- 
tial is included in consideration of the variance. 

A perspective on nonlinear I-V property. — Let us leave 
Z(x,t) with the Active field x and time t and return to 
y (n, N) with n intersections of neighboring segmented TDs 
and array length N. In the long array limit, the distribution 
converges to a Gaussian with variance N/12. 

Finally, we note the connection of TD to the nonlinearity in 
the I-V behavior. One can describe the average I-V property 
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/(V) := I(V, {q}), where the overline indicates the average 
for all sets {q}. In general, the offset charge distribution af- 
fects not only the value of the threshold, but also the trajectory 
of the electron between positive and negative electrodes. Each 
I(V, {qi}) is linear just above its threshold ||6|] as 

I{V, {q}) = G({q})(V - V*{{q}))H(V - V*({q})), (13) 

where %{-) denotes the Heaviside step function. The coef- 
ficient G depends on the trajectory of an electron and con- 
sequently on {q}. Here, let us consider ID arrays, where G 
is regarded as a constant for all offset charge distributions; 
i.e., the offset charge distribution influences only the value of 
the threshold. The average I-V property of ID arrays thus 
reduces to I lD (V) = J °° I(V, V th )P N (V th )dV th . Further, the 
conductance reduces to 

^=G £ P N (V th )dV th , (14) 

that is, the conductance is represented by the cumulative dis- 
tribution of Vth. 

In the model we employ, the conductance for long arrays is 
represented by the error function. Since it is not unusual that 
the TD is Gaussian, a conductance represented by the error 
function might be universal. In addition, in higher dimen- 
sional arrays, we can estimate an approximate I-V behavior 
by a superposition of ID paths, although it would be difficult 
to consider features such as meandering, bifurcation, and con- 
fluence. 

Summary. — We have obtained analytically the TD for a 
locally-coupled ID CB array containing N Coulomb islands. 
We first found an expression between V± and q\. Second, 
we introduced the segmented TD as a sum of products of the 
probability Ilv(fc) and the conditional probability Ujsr(n\k). 
Determining U]\r(n\k) leads to specific equations for the en- 
tire TD that perfectly describe our simulation results. In the 
long-array limit, the distribution converges to Gaussian form 
with variance N/12. In addition, we discussed a general char- 
acteristic of the nonlinear I-V behavior, where the cumulative 
distribution of the threshold voltage corresponds to the con- 
ductance. The current for each offset charge distribution and 
confirmation of this viewpoint will be discussed elsewhere. 

We also revealed that the distribution of the intersection is 
equivalent to the distribution Y(n,N), which indicates the 
probability for n — 1 upward steps for N — 1 aligned ob- 
jects of different height. Moreover, the distribution Z(x,t), 
which is equivalent to Y(n,N), satisfies a novel Fokker- 
Planck equation corresponding to active Brownian motion; 
i.e., overdamped motion of a Brownian particle in a harmonic 
potential that spreads with time. This relationship is analo- 
gous to Brownian motion and the binomial coefficients (i.e., 
the Pascal triangle). Further, the concept underlying the distri- 
bution of upward steps will be applicable to other nonequilib- 
rium and/or disordered systems. We focused on the derivation 
of the recurrence formula and the continuous limit in this Let- 



ter. It will be interesting to investigate characteristics of the 
novel Fokker-Planck equation. 
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